## libraries
library(DESeq2)
library(limma)
library(edgeR)
library(DT)
library(ComplexHeatmap)
#library(circlize)
library(corrplot)
library(survival)
#library(RColorBrewer)
# gene sets
data(human_gene_signatures)
ind_genes <- human_gene_signatures

# colors
data(color_palettes)

# variables
irf <- "Best Confirmed Overall Response"
ml <- "FMOne mutation burden per MB"
goi <- names(ind_genes)
tablesDir <- system.file("tables", 
  package="IMvigor210CoreBiologies")

# font sizes
labCex <- 0.9
namesCex <- 0.9
legendCex <- 0.9
titleCex <- 1
axisCex <- 0.9
titleF <- 1
# load
data(cds)  
cds2 <- cds
data(fmone)  
fmi <- fmone

# normalize 
geneNames <- setNames(fData(cds2)$Symbol, 
  as.character(rownames(fData(cds2))))
voomD <- filterNvoom(counts(cds2),
  minSamples=ncol(counts(cds2))/10,
  minCpm=0.25)
m <- voomD$E
m <- t(scale( t( m ),
  center=TRUE, 
  scale=TRUE)
)
# add signature scores to pData()
m2 <- m
rownames(m2) <- geneNames[rownames(m2)]

# calculate gene set scores
for (sig in goi) {
  pData(cds2)[, sig] <- NA
  genes <- ind_genes[[sig]]
  genes <- genes[genes %in% rownames(m2)]
  tmp <- m2[genes, , drop=FALSE]
  pData(cds2)[, sig] <- gsScore(tmp)
}
pData(cds2)$MKI67 <- m2["MKI67", ]

Figure 1a

oldMar <- par()$mar

par(mar=c(5.5, 4.1, 2, 0.5))

feat <- "IC Level"

tmpDat <- pData(cds2)[ !is.na(pData(cds2)[, feat]) & pData(cds2)[, irf] != "NE", ]
tmpDat[, irf] <- droplevels(tmpDat[, irf])
  
ic <- table(tmpDat[, feat], tmpDat[, "binaryResponse"])
pval <- signif(fisher.test(ic)$p.value, 2)
print(paste("Fisher P for IC by binary response:", pval))
## [1] "Fisher P for IC by binary response: 0.0038"
for (rr in c("CR", "PR")) {
  tmpDat$group <- ifelse(tmpDat[, irf] == rr, "R", "NR")
  tmpDat$group <- factor(tmpDat$group)
  ic <- table(tmpDat[, feat], tmpDat$group)
  print(rr)
  print("all")
  print(signif(fisher.test(ic)$p.value, 2))
  ic01 <- signif(fisher.test(ic[c("IC0", "IC1"),])$p.value, 2)
  ic02 <- signif(fisher.test(ic[c("IC0", "IC2+"),])$p.value, 2)
  ic12 <- signif(fisher.test(ic[c("IC1", "IC2+"),])$p.value, 2)
  ic2 <- rbind(ic01=colSums(ic[c("IC0", "IC1"),]),
    ic2=ic["IC2+",])
  ic2 <- signif(fisher.test(ic2)$p.value, 2)
  tmp <- c("IC0vsIC1"=ic01, "IC0vs2+"=ic02, "IC1vs2+"=ic12) 
  tmp <- ifelse(tmp * 3 > 1, 1, tmp * 3)
  print(c(tmp, "IC2vsOthers"=ic2))
}
## [1] "CR"
## [1] "all"
## [1] 0.0019
##    IC0vsIC1     IC0vs2+     IC1vs2+ IC2vsOthers 
##      1.0000      0.0126      0.0171      0.0006 
## [1] "PR"
## [1] "all"
## [1] 0.53
##    IC0vsIC1     IC0vs2+     IC1vs2+ IC2vsOthers 
##        1.00        0.93        1.00        0.30
ic <- table(tmpDat[, feat], tmpDat[, irf])
nSamples <- rowSums(ic)
ic <- prop.table(t(ic), 
  margin=2)
a <- barplot(ic,
  ylab="fraction of patients",
  cex.names=namesCex,
  cex.axis=axisCex,
  cex.lab=labCex,
  legend.text=rownames(ic),
  col=color_palettes$irf_palette,
  width=0.16,
  xlim=c(0,1),
  args.legend=list(bty="n",
     cex=legendCex,
     x="topright"),
  xaxt="n")
text(x = a, 
  y = par("usr")[3] - 0.03, 
  labels = levels(tmpDat[, feat]), 
  srt = -45, 
  xpd = TRUE,
  adj=0,
  cex=namesCex)
mtext(nSamples,
  side=3,
  at=a,
  line=0,
  cex=0.7)
mtext(paste("IC", "response", sep=", "),
  side=3,
  at=a[2],
  line=1,
  font=titleF,
  cex=titleCex)

par(mar=oldMar)

Figure 1b

addDots <- FALSE

orgMar <- par()$mar

par(mar=c(5.5, 4.1, 2, 0.5))

tmpDat <- pData(cds2)[pData(cds2)[, irf] != "NE",]
tmpDat[, irf] <- droplevels(tmpDat[, irf])

pval <- signif(wilcox.test(tmpDat[, "CD 8 T effector"] ~ tmpDat[, "binaryResponse"])$p.value, 
2)
print(paste("Wilcoxon test P for Teff by binary response:", pval))
## [1] "Wilcoxon test P for Teff by binary response: 0.0056"
pval <- signif(t.test(tmpDat[, "CD 8 T effector"] ~ tmpDat[, "binaryResponse"])$p.value, 
2)
print(paste("T test P for Teff by binary response:", pval))
## [1] "T test P for Teff by binary response: 0.0087"
yLab <- "CD8 T-effector score"
a <- boxplot(tmpDat[, "CD 8 T effector"] ~ tmpDat[, irf],
  ylab=yLab,
  col=color_palettes$irf_palette[levels(tmpDat[, irf])],
  cex.main=0.6,
  cex.axis=axisCex,
  cex.names=namesCex,
  cex.lab=labCex,
  ylim=c(-8, 6.5), 
  whisklty = 1)
if (addDots) {
  points(y=tmpDat[, "CD 8 T effector"],
    x=jitter(as.numeric(tmpDat[, irf]), factor=0.7),
    col=alpha("darkgrey", 0.7),
    pch=16)
}
mtext(a$n,
  side=3,
  at=1:nlevels(tmpDat[, irf]),
  line=0,
  cex=0.75)

mtext("CD8 T-effector, response",
  side=3,
  at=2.5,
  line=1,
  font=titleF,
  cex=titleCex)

ind <- tmpDat[, irf] %in% c("CR", "PR")
pvalCRpr <- t.test(tmpDat[ind, "CD 8 T effector"] ~ droplevels(tmpDat[ind, irf]))$p.value
ind <- tmpDat[, irf] %in% c("CR", "SD")
pvalCRsd <- t.test(tmpDat[ind, "CD 8 T effector"] ~ droplevels(tmpDat[ind, irf]))$p.value
ind <- tmpDat[, irf] %in% c("CR", "PD")
pvalCRpd <- t.test(tmpDat[ind, "CD 8 T effector"] ~ droplevels(tmpDat[ind, irf]))$p.value
pvalsCR <- c(CRpr=pvalCRpr, 
  CRsd=pvalCRsd, 
  CRpd=pvalCRpd)

ind <- tmpDat[, irf] %in% c("PR", "SD")
pvalPRsd <- t.test(tmpDat[ind, "CD 8 T effector"] ~ droplevels(tmpDat[ind, irf]))$p.value
ind <- tmpDat[, irf] %in% c("PR", "PD")
pvalPRpd <- t.test(tmpDat[ind, "CD 8 T effector"] ~ droplevels(tmpDat[ind, irf]))$p.value
ind <- tmpDat[, irf] %in% c("SD", "PD")
pvalSDpd <- t.test(tmpDat[ind, "CD 8 T effector"] ~ droplevels(tmpDat[ind, irf]))$p.value
pvalsOther <- c(PRsd=pvalPRsd, 
  PRsd=pvalPRpd, 
  SDpd=pvalSDpd)

pvalsAll <- c(pvalsCR, pvalsOther)
pvalsAll <- ifelse(pvalsAll * length(pvalsAll) > 1, 1, pvalsAll * length(pvalsAll))
print("t test P for all possible groupings:") 
## [1] "t test P for all possible groupings:"
print(pvalsAll)
##         CRpr         CRsd         CRpd         PRsd         PRsd         SDpd 
## 0.0386777173 0.0667884675 0.0002873751 1.0000000000 1.0000000000 0.2370005851
pvalsStar <- pLevel(pvalsAll)

segments(x0=1,
  x1=4,
  y0=6.45,
  y1=6.45,
  col="black",
  xpd=TRUE,
  lwd=1)
text(x=2.5,
  y=6.75,
  labels=pvalsStar["CRpd"],
  font=1,
  cex=1.2,
  xpd=TRUE)

segments(x0=1,
  x1=3,
  y0=6.06,
  y1=6.06,
  col="black",
  xpd=TRUE,
  lwd=1)
text(x=2,
  y=6.28,
  labels=pvalsStar["CRsd"],
  font=1,
  cex=1.6,
  xpd=TRUE)

segments(x0=1,
  x1=2,
  y0=5.5,
  y1=5.5,
  col="black",
  xpd=TRUE,
  lwd=1)
text(x=1.5,
  y=5.75,
  labels=pvalsStar["CRpr"],
  font=1,
  cex=1.2,
  xpd=TRUE)

par(mar=orgMar)

crAll <- factor(ifelse(tmpDat[, irf] == "CR", "CR", "nonCR"))
pvalCRall <- t.test(tmpDat[, "CD 8 T effector"] ~ crAll)$p.value
print(paste("t test P for CR vs all:", pvalCRall))
## [1] "t test P for CR vs all: 0.000205031608277682"
ind <- tmpDat[, irf] != "CR"
print("ANOVA P for PR/SD/PD")
## [1] "ANOVA P for PR/SD/PD"
getPfromAnova(tmpDat[ind, "CD 8 T effector"], droplevels(tmpDat[ind, irf]))
## [1] 0.0847

Figure 1c

bio <- "CD 8 T effector"
tmpDat <- pData(cds2)[!is.na(pData(cds2)[, bio]),]
tmpDat$group <- cut(tmpDat[, bio], 
  quantile(tmpDat[, bio], c(0, 0.25, 0.5, 0.75, 1)),
  include.lowest = TRUE)
tmpDat$group <- factor(tmpDat$group,
  labels=c("1st Quart", 
    "2nd Quart",
    "3rd Quart",
    "4th Quart"))
fittedSurv <- survfit(Surv(time=os, event=censOS) ~ group, 
              data=tmpDat)
diffSurv <- survdiff(Surv(time=os, event=censOS) ~ group,           
            data=tmpDat)
plotSurvival2(survFit=fittedSurv, 
    survDiff=diffSurv, 
    diff.factor=as.factor(tmpDat$group), 
    main="CD8 T-effector, survival", 
    cols=rev(color_palettes[["irf_palette"]][1:4]),
    ltypes=rep(1, nlevels(tmpDat$group)),
    pval="none",
    #coxphData=coxphDat,
    plotMedian=TRUE,
    plot.nrisk=FALSE,
    xLab="OS (months)",
    lwd=3,
    cexMedian=axisCex - 0.1,
    cexLegend=legendCex,
    cex.lab=labCex)

print(coxph(Surv(os, censOS) ~ 
  as.integer(tmpDat$group),
  data=tmpDat))
## Call:
## coxph(formula = Surv(os, censOS) ~ as.integer(tmpDat$group), 
##     data = tmpDat)
## 
##                              coef exp(coef) se(coef)      z       p
## as.integer(tmpDat$group) -0.15123   0.85965  0.05803 -2.606 0.00916
## 
## Likelihood ratio test=6.81  on 1 df, p=0.009064
## n= 348, number of events= 232

Figure 1d

addDots <- FALSE
alternative <- TRUE

## response by Mutation burden

oldMar <- par()$mar

par(mar=c(5.5, 4.1, 2, 0.5))

bio <- ml
feat <- irf

clinDat2 <- pData(cds2)[!is.na(pData(cds2)[, bio]) & pData(cds2)[, bio] > 0,]

tmpDat <- clinDat2[!is.na(clinDat2[, feat]),]
tmpDat <- tmpDat[tmpDat[, feat] != "NE",]
tmpDat[, feat] <- droplevels(tmpDat[, feat])
ind <- tmpDat[, bio] > 0
tmpDat[, bio] <- log2(tmpDat[, bio])

yLab <- "TMB"

a <- boxplot(tmpDat[ind, bio] ~ tmpDat[ind, feat],
  ylab=yLab,
  col=color_palettes$irf_palette,
  cex.main=0.6,
  cex.axis=axisCex,
  cex.names=namesCex,
  cex.lab=labCex,  
  whisklty = 1,
  ylim=c(1,log2(80)),
  yaxt="n")
axis(2,
  at=log2(c(1,2,5,10,20,50)),
  labels=c(1,2,5,10,20,50),
  cex=axisCex)
if (addDots) {
  points(y=tmpDat[ind, bio],
    x=jitter(as.numeric(tmpDat[ind, feat]), factor=0.7),
    col=alpha("darkgrey", 0.7),
    pch=16)
} 
  mtext(a$n,
  side=3,
  at=1:nlevels(tmpDat[, feat]),
  line=0,
  cex=0.8)

pval <- signif(wilcox.test(as.numeric(tmpDat[ind, bio]) ~ tmpDat[ind, "binaryResponse"])$p.value, 
    2)
print(paste("Wilcoxon Test p TMB by response:", pval))
## [1] "Wilcoxon Test p TMB by response: 1.4e-07"
pval <- signif(t.test(as.numeric(tmpDat[ind, bio]) ~ tmpDat[ind, "binaryResponse"])$p.value, 
    2)
print(paste("T Test p TMB by response:", pval))
## [1] "T Test p TMB by response: 6.9e-07"
mtext("TMB, response",
  side=3,
  at=2.5,
  line=1,
  font=titleF,
  cex=titleCex)

if (alternative) {
    yrange <- par("usr")[4] - par("usr")[3]
    yunit <- yrange/60
    segments(x0=1,
      x1=2,
      y0=par("usr")[4] - yunit * 4,
      y1=par("usr")[4] - yunit * 4,
      col="black",
      xpd=TRUE,
      lwd=1)
    segments(x0=3,
      x1=4,
      y0=par("usr")[4] - yunit * 4,
      y1=par("usr")[4] - yunit * 4,
      col="black",
      xpd=TRUE,
      lwd=1)
    segments(x0=1.5,
      x1=1.5,
      y0=par("usr")[4] - yunit * 4,
      y1=par("usr")[4] - yunit * 2,
      col="black",
      xpd=TRUE,
      lwd=1)
    segments(x0=3.5,
      x1=3.5,
      y0=par("usr")[4] - yunit * 4,
      y1=par("usr")[4] - yunit * 2,
      col="black",
      xpd=TRUE,
      lwd=1)
    segments(x0=1.5,
      x1=3.5,
      y0=par("usr")[4] - yunit * 2,
      y1=par("usr")[4] - yunit * 2,
      col="black",
      xpd=TRUE,
      lwd=1)
    text(y=par("usr")[4] - yunit ,
      x=2.5,
      labels=pLevel(pval),
      cex=axisCex) 
}

par(mar=oldMar)

Figure 1e

tmpDat <- pData(cds2)[!is.na(pData(cds2)[, ml]),]
tmpDat$group <- cut(tmpDat[, ml], 
  quantile(tmpDat[, ml], c(0, 0.25, 0.5, 0.75, 1)),
  include.lowest = TRUE)
tmpDat$group <- factor(tmpDat$group,
  labels=c("1st Quart", 
    "2nd Quart",
    "3rd Quart",
    "4th Quart"))
fittedSurv <- survfit(Surv(time=os, event=censOS) ~ group, 
  data=tmpDat)
diffSurv <- survdiff(Surv(time=os, event=censOS) ~ group,           
  data=tmpDat)
plotSurvival2(survFit=fittedSurv, 
  survDiff=diffSurv, 
  diff.factor=as.factor(tmpDat$group), 
  main="TMB, survival", 
  cols=rev(color_palettes[["irf_palette"]][1:4]),
  ltypes=rep(1, nlevels(tmpDat$group)),
  pval="none",
  #coxphData=coxphDat,
  plotMedian=TRUE,
  plot.nrisk=FALSE,
  xLab="OS (months)",
  lwd=3,
  cexMedian=axisCex,
  cexLegend=legendCex,
  cex.lab=labCex)

print(coxph(Surv(os, censOS) ~ 
  as.integer(tmpDat$group),
  data=tmpDat))
## Call:
## coxph(formula = Surv(os, censOS) ~ as.integer(tmpDat$group), 
##     data = tmpDat)
## 
##                              coef exp(coef) se(coef)      z     p
## as.integer(tmpDat$group) -0.29531   0.74430  0.06924 -4.265 2e-05
## 
## Likelihood ratio test=19  on 1 df, p=1.309e-05
## n= 272, number of events= 173

Figure 1f

ind <- !is.na(pData(fmi)$binaryResponse)
pData(fmi)$binaryResponse <- factor(pData(fmi)$binaryResponse,
  levels=c("CR/PR", "SD/PD"))

par(mar=c(5.5, 4.1, 2, 0))

path <- "DDR"

ids <- ind_genes[[path]]
ids <- ids[!ids %in% "TP53"]
tmp <- fmi[featureNames(fmi) %in% ids, ]
mutStatus <- any_mutation(tmp)
mutStatus <- colSums(mutStatus)

pval <- signif(fisher.test(table(pData(fmi)[ind, "binaryResponse"], mutStatus[ind] > 0))$p.value)
print(paste("Fisher P binary response by DDR mutation status:", pval))
## [1] "Fisher P binary response by DDR mutation status: 0.0116883"
b <- table(responder = fmi$binaryResponse[ind],
        mutant = factor(mutStatus[ind] > 0))
nSamples <- colSums(b)
d <- prop.table(b, 
margin=2)

a <- barplot(d,
  legend.text=colnames(d),
  col=color_palettes$response_palette[rownames(d)],
  cex.names=namesCex,
  cex.axis=axisCex,
  cex.lab=labCex,
  width=0.2,
  xlim=c(0,1),
  ylab="fraction of patients",
  xaxt="n",
  args.legend=list(x="topright",
    bty="n",
    fill=color_palettes$response_palette,
    legend=c("CR/PR", "SD/PD"),
    cex=legendCex))
text(x = a, 
  y = par("usr")[3] - 0.06, 
  labels = ifelse(colnames(d) == FALSE, "non-mutant", "mutant"), 
  srt = -45, 
  xpd = TRUE,
  adj=0,
  cex=namesCex,
  xpd=TRUE)
mtext(nSamples,
  side=3,
  at=a,
  line=0,
  cex=0.7)
mtext(paste(path, "w/o TP53, response"),
  side=3,
  at=mean(a),
  line=1,
  font=titleF,
  cex=titleCex)

mtext(pLevel(pval),
  at=mean(a),
  side=1,
  line=0.05,
  cex=namesCex,
  xpd=TRUE)  
segments(x0=a[1] - 0.05,
  x1=a[2] + 0.05,
  y0=-0.02,
  y1=-0.02,
  col="black",
  xpd=TRUE,
  lwd=1)

par(mar=oldMar)

Figure 1g

Pathway stats

ind <- !is.na(pData(fmi)$binaryResponse)
pathFMIout <- lapply(c("DDR", "CC Reg"), function(path) {
  ids <- ind_genes[[path]]
  tmp <- fmi[featureNames(fmi) %in% ids, ]
  mutStatus <- any_mutation(tmp)
  mutStatus <- colSums(mutStatus)
  t <- table(responder = fmi$binaryResponse[ind],
            mutant = factor(mutStatus[ind] > 0))
  f <- fisher.test( t )
  setNames(c( path, sum(mutStatus[ind] > 0), "response", f$estimate, f$p.value ),
    c("Pathway", "n mutant", "category", "estimate", "PVal"))
})
pathFMIout <- do.call(rbind, pathFMIout)
pathFMIout <- as.data.frame(pathFMIout,
  stringsAsFactors=FALSE)
pathFMIout$estimate <- signif(as.numeric(pathFMIout$estimate))
pathFMIout$PVal <- signif(as.numeric(pathFMIout$PVal))
#pathFMIout$AdjP <- p.adjust(pathFMIout$PVal, method="BH")

ind <- !is.na(pData(fmi)[,ml])
pathFMIout2 <- lapply(c("DDR", "CC Reg"), function(path) {
  ids <- ind_genes[[path]]
  tmp <- fmi[featureNames(fmi) %in% ids, ]
  mutStatus <- any_mutation(tmp)
  mutStatus <- colSums(mutStatus)
  f <- t.test( as.numeric(pData(fmi)[ind,ml]) ~ mutStatus[ind] > 0, 
    conf.int = TRUE )
  setNames(c( path, sum(mutStatus[ind] > 0), "TMB", f$estimate[1] - f$estimate[2], f$p.value ),
    c("Pathway", "n mutant", "category", "estimate", "PVal"))
})
pathFMIout2 <- do.call(rbind, pathFMIout2)
pathFMIout2 <- as.data.frame(pathFMIout2,
  stringsAsFactors=FALSE)
pathFMIout2$estimate <- signif(as.numeric(pathFMIout2$estimate))
pathFMIout2$PVal <- signif(as.numeric(pathFMIout2$PVal))
#pathFMIout2$AdjP <- p.adjust(pathFMIout2$PVal, method="BH")
pathFMIout <- rbind(pathFMIout, pathFMIout2)
datatable(pathFMIout,
  rownames=FALSE)

Single gene stats

tmp <- any_mutation(fmi)
tmp <- tmp[rowSums(tmp) > 0, ]
pathFMIout3 <- lapply(rownames(tmp), function(gene) {
  mutStatus <- tmp[gene,]
  ind <- !is.na(pData(fmi)[,ml])
  f <- try(t.test( as.numeric(pData(fmi)[ind,ml]) ~ mutStatus[ind] > 0, 
    conf.int = TRUE ))
  if (!is(f, "try-error")) {
    a <- setNames(c( gene, 
        sum(mutStatus[ind] > 0), 
        "TMB", 
        f$estimate[1]-f$estimate[2], 
        f$p.value ),
      c("Gene", "n mutant", "category", "estimate", "PVal"))
  } else {
    a <- setNames(c( gene, 
        sum(mutStatus[ind] > 0), 
        "TMB", 
        NA, 
        NA ),
      c("Gene", "n mutant", "category", "estimate", "PVal"))
  }
  ind <- !is.na(pData(fmi)$binaryResponse)
  t <- table(responder = fmi$binaryResponse[ind],
            mutant = factor(mutStatus[ind] > 0, levels=c("TRUE", "FALSE")))
  f <- fisher.test( t )
  b <- setNames(c( gene, sum(mutStatus[ind] > 0), "response", f$estimate, f$p.value ),
    c("Gene", "n mutant", "category", "estimate", "PVal"))
  rbind(a, b)
})

Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations Error in t.test.default(x = c(2.7, 6.31, 5.41, 9.91, 9.01, 11.71, 14.41, : not enough ‘y’ observations

pathFMIout3 <- do.call(rbind, pathFMIout3)
pathFMIout3 <- as.data.frame(pathFMIout3,
  stringsAsFactors=FALSE)
pathFMIout3$estimate <- signif(as.numeric(pathFMIout3$estimate))
pathFMIout3$PVal <- signif(as.numeric(pathFMIout3$PVal))

## adjust p-values
pathFMIout3$AdjP <- NA
pathFMIout3$AdjP[pathFMIout3$category == "response"] <- p.adjust(
  pathFMIout3$PVal[pathFMIout3$category == "response"] , method="BH")
pathFMIout3$AdjP[pathFMIout3$category == "TMB"] <- p.adjust(
  pathFMIout3$PVal[pathFMIout3$category == "TMB"] , method="BH")
pathFMIout3$inDDR <- pathFMIout3$Gene %in% ind_genes$DDR
pathFMIout3$inCCReg <- pathFMIout3$Gene %in% ind_genes$"CC Reg"
datatable(pathFMIout3,
  rownames=FALSE)
path <- "DDR"

fmiPathP <- pathFMIout[pathFMIout$category == "TMB", ]
fmiGeneP <- pathFMIout3[pathFMIout3$category == "TMB", ]

ids <- ind_genes[[path]]
ids <- ids[ids %in% featureNames(fmi)]
fmiGoi <- fmi[featureNames(fmi) %in% ids, ]
sOrder <- order(as.numeric(pData(fmiGoi)[, ml]), decreasing=TRUE)

mutGoi <- any_mutation(fmiGoi)

tmp1 <- mutGoi[rowSums(mutGoi) > 0,]
tmp2 <- colSums(tmp1)
tmp2 <- tmp2 > 0
if ("TP53" %in% ids) {
  tmp3 <- mutGoi[ids[! ids %in% "TP53"], ]
  tmp3 <- colSums(tmp3)
  tmp3 <- tmp3 > 0
  tmp2 <- rbind(tmp2, tmp3)
}
toSort <- rowSums(tmp1)
toSort <- sort(toSort, decreasing=TRUE)
mutGoi <- rbind(tmp1[names(toSort),], tmp2)
rownames(mutGoi)[rownames(mutGoi) == "tmp2"] <- "pathway"
if ("TP53" %in% ids) {
  rownames(mutGoi)[rownames(mutGoi) == "tmp3"] <- "pathway w/o TP53"
}


mutGoi <- apply(mutGoi, 2, function(x) {
  ifelse(x, "mutant", "")
})

fmiCols <- c(mutant="black")
alter_fun = function(x, y, w, h, v) {
    if(v["mutant"]) grid.rect(x, y, w*0.9, h*0.9, gp = gpar(fill = fmiCols["mutant"], 
      col = NA))
  }


ha = HeatmapAnnotation(TMB = anno_barplot(as.numeric(pData(fmiGoi)[sOrder, ml]),
    border=FALSE,
    gp = gpar(fill="black")),
  annotation_legend_param=list(gp = gpar(fontsize = 1)),
  show_annotation_name = TRUE,
  annotation_name_side="left",
  annotation_name_gp = gpar(fontsize = 10))

rownames(mutGoi) <- paste0(rownames(mutGoi),
  ifelse(rownames(mutGoi) %in% fmiGeneP$Gene[fmiGeneP$AdjP < 0.05], "*", ""))

op1 <- oncoPrint(mutGoi[, sOrder], 
  col=fmiCols,
  alter_fun=alter_fun,
  row_order=NULL,
  column_order=NULL,
  column_title = "Mutations in cell cycle regulation genes", 
  top_annotation=ha,
  show_pct=TRUE,
  pct_gp = gpar(fontsize = 9),
  row_names_gp = gpar(fontsize = 9),
  column_title_gp = gpar(fontsize = 10),
  heatmap_legend_param=list(title ="Mutation Status"),
  show_row_barplot = FALSE,
  show_heatmap_legend = FALSE)


op2 <- oncoPrint(mutGoi[, sOrder], 
  col=fmiCols,
  alter_fun=alter_fun,
  row_order=NULL,
  column_order=NULL,
  column_title = "Mutations in DDR genes", 
  top_annotation=ha,
  show_pct=TRUE,
  pct_gp = gpar(fontsize = 9),
  row_names_gp = gpar(fontsize = 9),
  column_title_gp = gpar(fontsize = 10),
  heatmap_legend_param=list(title ="Mutation Status"),
  show_row_barplot = FALSE,
  show_heatmap_legend = FALSE)

draw(op2, 
  padding = unit(c(1, 5, 0, 1), "mm"),
  row_title="mutation prevalence",
  row_title_gp=gpar(fontsize = 10))

Figure 1h

plotDots <- FALSE
binary <- FALSE

orgMar <- par()$mar

tmpDat <- pData(cds2)
tmpDat$TGFB1 <- scale(m2["TGFB1", ], center=TRUE, scale=TRUE)

tmpDat <- tmpDat[tmpDat[, irf] != "NE",]
tmpDat[, irf] <- droplevels(tmpDat[, irf])

sig <- "TGFB1"
  
par(mar=c(5.5, 4.1, 2, 0.5))

pval <- signif(wilcox.test(tmpDat[, sig] ~ tmpDat[, "binaryResponse"])$p.value, 
  2)
print(paste("Wilcoxon test P", sig, "score by binary Response:", pval))
## [1] "Wilcoxon test P TGFB1 score by binary Response: 5.6e-05"
pval <- signif(t.test(tmpDat[, sig] ~ tmpDat[, "binaryResponse"])$p.value, 
  2)
print(paste("T test P", sig, "score by binary Response:", pval))
## [1] "T test P TGFB1 score by binary Response: 0.00011"
yLab <- paste(sig, "expression")

a <- boxplot(tmpDat[, sig] ~ tmpDat[, irf],
  ylab=yLab,
  col=color_palettes[["irf_palette"]][levels(tmpDat[, irf])],
  cex.axis=axisCex,
  cex.names=namesCex,
  cex.lab=labCex, 
  ylim=c(-5, 5),
  whisklty = 1)
if (plotDots) {
  points(y=tmpDat[, sig],
    x=jitter(as.numeric(tmpDat[, irf]), factor=0.7),
    col=alpha("darkgrey", 0.7),
    pch=16)
}
mtext(a$n,
  side=3,
  at=1:nlevels(tmpDat[, irf]),
  line=0,
  cex=0.75)
mtext(paste0(sig, ", response"),
  side=3,
  at=2.5,
  line=1,
  font=titleF,
  cex=titleCex)

yrange <- par("usr")[4] - par("usr")[3]
yunit <- yrange/60
segments(x0=1,
  x1=2,
  y0=par("usr")[4] - yunit * 4,
  y1=par("usr")[4] - yunit * 4,
  col="black",
  xpd=TRUE,
  lwd=1)
segments(x0=3,
  x1=4,
  y0=par("usr")[4] - yunit * 4,
  y1=par("usr")[4] - yunit * 4,
  col="black",
  xpd=TRUE,
  lwd=1)
segments(x0=1.5,
  x1=1.5,
  y0=par("usr")[4] - yunit * 4,
  y1=par("usr")[4] - yunit * 2,
  col="black",
  xpd=TRUE,
  lwd=1)
segments(x0=3.5,
  x1=3.5,
  y0=par("usr")[4] - yunit * 4,
  y1=par("usr")[4] - yunit * 2,
  col="black",
  xpd=TRUE,
  lwd=1)
segments(x0=1.5,
  x1=3.5,
  y0=par("usr")[4] - yunit * 2,
  y1=par("usr")[4] - yunit * 2,
  col="black",
  xpd=TRUE,
  lwd=1)
text(y=par("usr")[4] - yunit ,
  x=2.5,
  labels=pLevel(pval),
  cex=ifelse(pLevel(pval) == "n.s.", axisCex - 0.1, axisCex))   

par(mar=orgMar)

Figure 1i

bio <- "TGFB1"

tmpDat  <- pData(cds2)
tmpDat$TGFB1 <- m2["TGFB1", ]

tmpDat <- tmpDat[!is.na(tmpDat[, bio]),]
tmpDat$group <- cut(tmpDat[, bio], 
  quantile(tmpDat[, bio], c(0, 0.25, 0.5, 0.75, 1)),
  include.lowest = TRUE)
tmpDat$group <- factor(tmpDat$group,
  labels=c("1st Quart", 
    "2nd Quart",
    "3rd Quart",
    "4th Quart"))
fittedSurv <- survfit(Surv(time=os, event=censOS) ~ group, 
              data=tmpDat)
diffSurv <- survdiff(Surv(time=os, event=censOS) ~ group,           
            data=tmpDat)
plotSurvival2(survFit=fittedSurv, 
    survDiff=diffSurv, 
    diff.factor=as.factor(tmpDat$group), 
    main="TGFB1, survival", 
    cols=rev(color_palettes[["irf_palette"]][1:4]),
    ltypes=rep(1, nlevels(tmpDat$group)),
    pval="none",
    #coxphData=coxphDat,
    plotMedian=TRUE,
    plot.nrisk=FALSE,
    xLab="OS (months)",
    lwd=3,
    cexMedian=axisCex - 0.1,
    cexLegend=legendCex,
    cex.lab=labCex)

print(coxph(Surv(os, censOS) ~ 
  as.integer(tmpDat$group),
  data=tmpDat))
## Call:
## coxph(formula = Surv(os, censOS) ~ as.integer(tmpDat$group), 
##     data = tmpDat)
## 
##                             coef exp(coef) se(coef)     z       p
## as.integer(tmpDat$group) 0.15298   1.16530  0.05905 2.591 0.00958
## 
## Likelihood ratio test=6.74  on 1 df, p=0.009433
## n= 348, number of events= 232

Extended Data

ED Figure 1

table(IC=pData(cds2)[, "IC Level"], TC=pData(cds2)[, "TC Level"])
##       TC
## IC     TC0 TC1 TC2+
##   IC0   93   4    0
##   IC1  111  10   11
##   IC2+  71   8   39
tmp <- pData(cds2)[!is.na(pData(cds2)[, "TC Level"]), ]
a <- ifelse(tmp[, "TC Level"] %in% "TC0", "TC-", "TC+")
b <- ifelse(tmp[, "IC Level"] %in% "IC0", "IC-", "IC+")
table(IC=b, TC=a)
##      TC
## IC    TC- TC+
##   IC-  93   4
##   IC+ 182  68
prop.table(table(IC=b, TC=a), 2)
##      TC
## IC           TC-        TC+
##   IC- 0.33818182 0.05555556
##   IC+ 0.66181818 0.94444444
par(mar=c(5.5, 4.1, 2, 0.5))

for (feat in "TC Level") {
  tmpDat <- pData(cds2)[ !is.na(pData(cds2)[, feat]) & pData(cds2)[, irf] != "NE", ]
  tmpDat[, irf] <- droplevels(tmpDat[, irf])
  ic <- table(tmpDat[, feat], tmpDat[, "binaryResponse"])
  pval <- signif(fisher.test(ic)$p.value, 2)
  print(paste("Fisher P for TC by binary response:", pval))
  ic <- table(tmpDat[, feat], tmpDat[, irf])
  nSamples <- rowSums(ic)
  ic <- prop.table(t(ic), 
    margin=2)
  a <- barplot(ic,
    #main=paste("P:", pval),
    ylab="fraction of patients",
    #xlab=feat,
    cex.names=namesCex,
    cex.axis=axisCex,
    cex.lab=labCex,
    legend.text=rownames(ic),
    col=color_palettes[["irf_palette"]],
    width=0.16,
    xlim=c(0,1),
    #las=2,
    args.legend=list(bty="n",
          #bg=alpha("white", 0.7),
          #box.col=alpha("white", 0.7),
          cex=1.1,
          x="topright"),
    xaxt="n")
  text(x = a, 
    y = par("usr")[3] - 0.03, 
    labels = levels(tmpDat[, feat]), 
    srt = -45, 
    #pos = 1, 
    xpd = TRUE,
    adj=0,
    cex=namesCex)
  mtext(nSamples,
    side=3,
    at=a,
    line=0,
    cex=0.7)
  #mtext(paste("P:", pval),
  #  side=1,
  #  at=a[2],
  #  line=3.5,
  #  cex=0.8)
  mtext(paste("TC", "response", sep=", "),
    side=3,
    at=a[2],
    line=1,
    font=titleF,
    cex=titleCex)
}
## [1] "Fisher P for TC by binary response: 0.72"

oldMar <- par()$mar

par(mar=c(5.5, 3.5, 2, 0.5))

icAss <- read.csv(file.path(tablesDir,
    "SupplementaryTableS10.csv"),
  stringsAsFactor=FALSE,
  check.names=FALSE)
icAss <- icAss[, 1:8]
icAss <- icAss[!is.na(icAss$Symbol), ]
interferome <- read.table(file.path(tablesDir,
    "20170503081937_GeneSearchResults_InterferomeTop300IC.txt"),
  stringsAsFactor=FALSE,
  sep="\t",
  check.names=FALSE,
  skip=19,
  header=TRUE)

cols <- "darkgrey"

plot(x=icAss[,"logFC"],
  y=-log10(icAss[,"P (adj.)"]),
  pch=16,
  cex=0.8,
  col=cols,
  ylab="",
  xlab="",
  lwd=2,
  cex.axis=axisCex,
  cex.lab=labCex)
mtext(expression(paste("-log"["10"], " adj. p")),
  side=2,
  line=2,
  cex=axisCex)
abline(h=-log10(0.05), lty=2, col="black")
gs <- list("Interferome"=interferome$"Gene Name"[interferome$"Gene Name" %in% icAss$Symbol[1:307]],
  "CD 8 T effector"=ind_genes[["CD 8 T effector"]],
  "Immune Checkpoint"=ind_genes[["Immune Checkpoint"]])
gsC <- list("Interferome"="orange1",
  "CD 8 T effector"="magenta2",
  "Immune Checkpoint"="olivedrab3")
for (g in names(gs)) {
  ind <- which(icAss$Symbol %in% gs[[g]])
  points(x=icAss[ind,"logFC"],
    y=-log10(icAss[ind,"P (adj.)"]),
    pch=16,
    cex=0.8,
    col=gsC[[g]])
}
legend("topleft",
  pch=16,
  col=c("orange1", "magenta2", "olivedrab3"),
  legend=c("IFNg inducible", "CD8 T-effector", "Immune checkpoint"),
  bty = "n",
  cex=0.6)
labelSub <- icAss[icAss$Symbol %in% 
  c(ind_genes[["CD 8 T effector"]], ind_genes[["Immune Checkpoint"]]), ]
labelPos <- c(2, 
  2,
  2, 
  4, 
  rep(3, 2), 
  4,
  3,
  4,
  2,
  1,
  2,
  1,
  2, 
  1)
text(x=labelSub[,"logFC"], 
  y=-log10(labelSub[,"P (adj.)"]), 
  labels=labelSub$Symbol, 
  cex= 0.4, 
  pos=labelPos)
mtext("Expression, IC",
    side=3,
    at=0.5,
    line=1,
    font=titleF,
    cex=titleCex)
mtext("log fold change",
    side=1,
    at=0.5,
    line=2,
    cex=axisCex)

par(mar=oldMar)
addDots <- FALSE

sig <- "CD 8 T effector"
feat <- "IC Level"

oldMar <- par()$mar

par(mar=c(5.5, 4.1, 2, 2))

tmpDat <- pData(cds2)[ !is.na(pData(cds2)[, feat]), ]
a <- boxplot(tmpDat[, sig] ~ tmpDat[, feat],
  ylab="CD8 T-effector score",
  col=color_palettes$ic_palette[levels(tmpDat[, feat])],
  cex.axis=axisCex,
  cex.lab=labCex,
  cex.names=namesCex, 
  whisklty = 1,
  xaxt="n")
axis(1,
  at=1:nlevels(tmpDat[, feat]),
  labels=rep("", nlevels(tmpDat[, feat])))
text(x = 1:nlevels(tmpDat[, feat]), 
  y = par("usr")[3] - 1.5, 
  labels = levels(tmpDat[, feat]), 
  xpd = TRUE,
  cex=namesCex)

if (addDots) {
  points(y=tmpDat[, sig],
    x=jitter(as.numeric(tmpDat[, feat]), factor=0.7),
    col=alpha("darkgrey", 0.7),
    pch=16)
}

mtext(a$n,
  side=3,
  at=1:nlevels(tmpDat[, feat]),
  line=0,
  cex=0.8)
mtext(paste("CD8 T-effector", "IC", sep=", "),
  side=3,
  at=2,
  line=1,
  font=titleF,
  cex=titleCex)

pval <- signif(getPfromAnova(tmpDat[, sig], tmpDat[, feat]), 2)
print(paste("ANOVA P for Teff signature by IC level:", pval))
## [1] "ANOVA P for Teff signature by IC level: 4.2e-35"
par(mar=oldMar)
addDots <- FALSE
alternative <- TRUE

## response by Mutation burden

oldMar <- par()$mar

par(mar=c(5.5, 4.1, 2, 0.5))

bio <- "Neoantigen burden per MB"
feat <- irf

clinDat2 <- pData(cds2)[!is.na(pData(cds2)[, bio]) & pData(cds2)[, bio] > 0,]

tmpDat <- clinDat2[!is.na(clinDat2[, feat]),]
tmpDat <- tmpDat[tmpDat[, feat] != "NE",]
tmpDat[, feat] <- droplevels(tmpDat[, feat])
ind <- tmpDat[, bio] > 0
tmpDat[, bio] <- log2(tmpDat[, bio])

yLab <- "TNB"

a <- boxplot(tmpDat[ind, bio] ~ tmpDat[ind, feat],
  ylab=yLab,
  col=color_palettes$irf_palette,
  cex.main=0.6,
  cex.axis=axisCex,
  cex.names=namesCex,
  cex.lab=labCex,  
  whisklty = 1,
  ylim=c(-5, log2(30)),
  yaxt="n")
axis(2,
  at=log2(c(0.05,0.20,1.00,5.00,20)),
  labels=round(c(0.05,0.20,1.00,5.00,20.00), 2),
  cex=axisCex)
if (addDots) {
  points(y=tmpDat[ind, bio],
    x=jitter(as.numeric(tmpDat[ind, feat]), factor=0.7),
    col=alpha("darkgrey", 0.7),
    pch=16)
} 
  mtext(a$n,
  side=3,
  at=1:nlevels(tmpDat[, feat]),
  line=0,
  cex=0.8)

pval <- signif(wilcox.test(as.numeric(tmpDat[ind, bio]) ~ tmpDat[ind, "binaryResponse"])$p.value, 
    2)
print(paste("Wilcoxon p-value, neoantigens by response:", pval))
## [1] "Wilcoxon p-value, neoantigens by response: 5.3e-09"
pval <- signif(t.test(as.numeric(tmpDat[ind, bio]) ~ tmpDat[ind, "binaryResponse"])$p.value, 
    2)
print(paste("T test p-value, neoantigens by response:", pval))
## [1] "T test p-value, neoantigens by response: 2.7e-09"
mtext("Neoantigens, response",
  side=3,
  at=2.5,
  line=1,
  font=titleF,
  cex=titleCex)

if (alternative) {
    yrange <- par("usr")[4] - par("usr")[3]
    yunit <- yrange/60
    segments(x0=1,
      x1=2,
      y0=par("usr")[4] - yunit * 4,
      y1=par("usr")[4] - yunit * 4,
      col="black",
      xpd=TRUE,
      lwd=1)
    segments(x0=3,
      x1=4,
      y0=par("usr")[4] - yunit * 4,
      y1=par("usr")[4] - yunit * 4,
      col="black",
      xpd=TRUE,
      lwd=1)
    segments(x0=1.5,
      x1=1.5,
      y0=par("usr")[4] - yunit * 4,
      y1=par("usr")[4] - yunit * 2,
      col="black",
      xpd=TRUE,
      lwd=1)
    segments(x0=3.5,
      x1=3.5,
      y0=par("usr")[4] - yunit * 4,
      y1=par("usr")[4] - yunit * 2,
      col="black",
      xpd=TRUE,
      lwd=1)
    segments(x0=1.5,
      x1=3.5,
      y0=par("usr")[4] - yunit * 2,
      y1=par("usr")[4] - yunit * 2,
      col="black",
      xpd=TRUE,
      lwd=1)
    text(y=par("usr")[4] - yunit ,
      x=2.5,
      labels=pLevel(pval),
      cex=axisCex) 
}

par(mar=oldMar)
tmpDat <- pData(cds2)[!is.na(pData(cds2)[, "Neoantigen burden per MB"]),]
tmpDat$group <- cut(tmpDat[, "Neoantigen burden per MB"], 
  quantile(tmpDat[, "Neoantigen burden per MB"], c(0, 0.25, 0.5, 0.75, 1)),
  include.lowest = TRUE)
tmpDat$group <- factor(tmpDat$group,
  labels=c("1st Quart", 
    "2nd Quart",
    "3rd Quart",
    "4th Quart"))
fittedSurv <- survfit(Surv(time=os, event=censOS) ~ group, 
  data=tmpDat)
diffSurv <- survdiff(Surv(time=os, event=censOS) ~ group,           
  data=tmpDat)
plotSurvival2(survFit=fittedSurv, 
  survDiff=diffSurv, 
  diff.factor=as.factor(tmpDat$group), 
  main="TNB, survival", 
  cols=rev(color_palettes[["irf_palette"]][1:4]),
  ltypes=rep(1, nlevels(tmpDat$group)),
  pval="none",
  #coxphData=coxphDat,
  plotMedian=TRUE,
  plot.nrisk=FALSE,
  xLab="OS (months)",
  lwd=3,
  cexMedian=axisCex - 0.1,
  cexLegend=legendCex,
  cex.lab=labCex)

print(coxph(Surv(os, censOS) ~ 
  as.integer(tmpDat$group),
  data=tmpDat))
## Call:
## coxph(formula = Surv(os, censOS) ~ as.integer(tmpDat$group), 
##     data = tmpDat)
## 
##                              coef exp(coef) se(coef)      z        p
## as.integer(tmpDat$group) -0.26824   0.76472  0.07072 -3.793 0.000149
## 
## Likelihood ratio test=14.75  on 1 df, p=0.0001229
## n= 245, number of events= 160
print(coxph(Surv(os, censOS) ~ 
  as.ordered(tmpDat$group),
  data=tmpDat))
## Call:
## coxph(formula = Surv(os, censOS) ~ as.ordered(tmpDat$group), 
##     data = tmpDat)
## 
##                               coef exp(coef) se(coef)      z        p
## as.ordered(tmpDat$group).L -0.6645    0.5145   0.1735 -3.829 0.000129
## as.ordered(tmpDat$group).Q -0.3235    0.7236   0.1649 -1.962 0.049800
## as.ordered(tmpDat$group).C -0.2942    0.7451   0.1566 -1.878 0.060337
## 
## Likelihood ratio test=21.28  on 3 df, p=9.203e-05
## n= 245, number of events= 160
keggRes <- read.csv(file.path(tablesDir,
    "SupplementaryTableS1.csv"),
  stringsAsFactor=FALSE,
  check.names=FALSE)
keggRes <- keggRes[, 1:6]
keggRes <- keggRes[!is.na(keggRes$S), ]
ind <- rev(1:nrow(keggRes))
keggRes <- keggRes[ind,]

up <- keggRes$Direction == "Up"
keggRes$"P (adj.)"[up] <- log10(keggRes$"P (adj.)"[up]) * -1
down <- keggRes$Direction == "Down"
keggRes$"P (adj.)"[down] <- log10(keggRes$"P (adj.)"[down])

resC1 <- read.csv(file.path(tablesDir,
    "SupplementaryTableS2.csv"),
  stringsAsFactor=FALSE,
  check.names=FALSE)

keggRes$RepGenes <- NA_character_
for (i in 1:nrow(keggRes)) {
  intGenes <- unlist(strsplit(keggRes$"Identified Genes"[i], ", "))
  if (length(intGenes) > 7) {
    intGenes <- resC1$Symbol[resC1$Symbol %in% intGenes][1:7]
    intGenes <- c(intGenes, "...")
  }
  intGenes <- paste(intGenes, collapse=", ")
  keggRes$RepGenes[i] <- intGenes
}

layout(matrix(c(1,1,0,0,0,0,1,1,0,0,0,0), 2, 6, byrow=TRUE))

proliferationPathways <- c("DNA replication", 
  "Cell cycle", 
  "Systemic lupus erythematosus", 
  "Alcoholism",
  "Oocyte meiosis",
  "Viral carcinogenesis",
  "Progesterone-mediated oocyte maturation")
ddrPathways <- c("Fanconi anemia pathway", "Homologous recombination", "Mismatch repair", "Nucleotide excision repair", "Base excision repair", "p53 signaling pathway")

setCols <- ifelse(keggRes$Name %in% proliferationPathways,
  "turquoise4",
    ifelse(keggRes$Name %in% ddrPathways,
      "magenta4", 
        ifelse(keggRes$Name %in% "Pathways in cancer", "orange2", "darkgrey")))

a <- barplot(keggRes$"P (adj.)",
  horiz=TRUE,
  xlab="",
  xlim=c(-2, 6),
  col=setCols,
  axes=FALSE,
  cex.main=0.9)
mtext("Kegg pathway",
  side=4,
  at=max(a[,1]) + 1.5,
  las=2,
  line=3,
  cex=titleCex)
mtext("Top-scoring genes",
  side=4,
  at=max(a[,1]) + 1.5,
  las=2,
  line=24,
  cex=titleCex)
mtext(side=4,
  at=a[,1],
  text=paste0(keggRes$Name,
    ifelse(keggRes$rep, "*", "")),
  las=2,
  line=1.5,
  cex=(namesCex - 0.1))
mtext(keggRes$RepGenes,
  side=4,
  at=a[,1],
  line=22,
  cex=(namesCex - 0.1),
  las=2)
axis(side=1,
  at=c(-2, seq(from=2, to=6, by=2)),
  labels=c(2, seq(from=2, to=6, by=2)),
  tick=TRUE,
  xpd=TRUE,
  cex.axis=0.9)
segments(x0=0,
  x1=0,
  y0=-0.6,
  y1=max(a[,1])+0.7,
  col="black",
  xpd=TRUE,
  lwd=1)
mtext(expression(paste("-log"["10"], " adj. p")),
  at=2,
  side=1,
  line=2.4,
  cex=0.9)

mtext(c("Positive association with TMB", "Negative"), 
  side=2, 
  line=1.2, 
  at=c(a[,1][9], mean(c(a[,1][1], a[,1][2])) + 0.2),
  cex=axisCex)
segments(x0=-2.1,
  x1=-2.1,
  y0=a[,1][4] - 0.3,
  y1=a[,1][length(a[,1])] + 0.3,
  col="black",
  xpd=TRUE,
  lwd=1)
segments(x0=-2.1,
  x1=-2.1,
  y0=a[,1][1] - 0.5,
  y1=a[,1][3] + 0.5,
  col="black",
  xpd=TRUE,
  lwd=1)

legend("bottomright",
  fill=c("turquoise4", "magenta4", "orange2"),
  legend=c("Proliferation", "DNA damage response", expression(paste("TGF-", beta))),
  bty="n",
  cex=legendCex + 0.3)

goi <- c("CD 8 T effector",
  "gene19",
  "Fanconi",
  "MKI67",
  "Cell cycle",
  "DNA replication",
  "Histones",
  "Homologous recombination",
  "Mismatch Repair",
  "Nucleotide excision repair",
  "CC Reg")

pLabels <- sub("CC Reg", "Cell cycle regulation", goi)
pLabels <- sub("CD 8 T e", "CD8 T-e", pLabels)

tmp <- pData(cds2)[, goi]
colnames(tmp) <- sub("CC Reg", "Cell cycle regulation", colnames(tmp))
colnames(tmp) <- sub("CD 8 T e", "CD8 T-e", colnames(tmp))
colnames(tmp) <- sub("CC Reg", "Cell cycle regulation", colnames(tmp))
colnames(tmp) <- sub("Fanconi", "Fanconi anemia pathway", colnames(tmp))
colnames(tmp) <- sub("p53", "p53-like", colnames(tmp))
colnames(tmp) <- sub("gene19", "Pan-F-TBRS", colnames(tmp))

M <- cor(tmp)
corrplot.mixed(M,
  lower="circle",
  upper="number",
  number.cex = 0.4,
  order="AOE",
  tl.cex = 0.9, 
  tl.col = "black",
  xpd=TRUE,
  tl.pos="lt",
  diag="n")

addDots <- FALSE

sig <- "APOBEC3A"

oldMar <- par()$mar


## plot APOBEC expression by response
tmpDat <- pData(cds2)
tmpDat[, sig] <- m2[sig, ]
tmpDat <- tmpDat[tmpDat[, irf] != "NE",]
tmpDat[, irf] <- droplevels(tmpDat[, irf])

pval <- signif(wilcox.test(tmpDat[, sig] ~ tmpDat[, "binaryResponse"])$p.value, 
2)
print(paste("Wilcoxon test P APOBEC3A expression by binary response:", pval))
## [1] "Wilcoxon test P APOBEC3A expression by binary response: 0.0071"
pval <- signif(t.test(tmpDat[, sig] ~ tmpDat[, "binaryResponse"])$p.value, 
2)
print(paste("T test P APOBEC3A expression by binary response:", pval))
## [1] "T test P APOBEC3A expression by binary response: 0.015"
yLab <- expression(paste("APOBEC3A expression (log"["2"], ")"))
par(mar=c(5.5, 4, 2, 0.5))

a <- boxplot(tmpDat[, sig] ~ tmpDat[, irf],
  ylab=yLab,
  col=color_palettes$irf_palette[levels(tmpDat[, irf])],
  cex.main=0.6,
  cex.axis=axisCex,
  cex.names=namesCex,
  cex.lab=labCex,  
  whisklty = 1)
if (addDots) {
  points(y=tmpDat[, sig],
    x=jitter(as.numeric(tmpDat[, irf]), factor=0.7),
    col=alpha("darkgrey", 0.7),
    pch=16)
}
mtext(a$n,
  side=3,
  at=1:nlevels(tmpDat[, irf]),
  line=0,
  cex=0.75)
mtext(paste(sig, "response", sep=", "),
  side=3,
  at=2.5,
  line=1,
  font=titleF,
  cex=titleCex)

## plot APOBEC expression by TMB
tmpDat <- pData(cds2)
tmpDat[, sig] <- m2[sig, ]
ind <- tmpDat[, ml] > 0
tmpDat[, ml] <- log2(tmpDat[, ml])

plot(tmpDat[ind, sig] ~ tmpDat[ind, ml],
    xlab="",
    ylab=expression(paste("APOBEC3A expression (log"["2"], ")")),
    pch=16,
    col="darkgrey",
    cex.axis=0.8,
    cex.lab=labCex,
    xaxt="n")
axis(1,
  at=log2(c(1,2,5,10,20,50)),
  labels=c(1,2,5,10,20,50),
  cex=axisCex)

fit <- lm(tmpDat[ind, sig] ~ tmpDat[ind, ml])
abline(fit, lty=2)

tmp <- cor.test(tmpDat[ind, sig], tmpDat[ind, ml], method="pearson")
text(paste("cor:", round(tmp$estimate, 2), 
    "\nP:", signif(tmp$p.value, 2)),
  font=2,
  cex=0.7,
  x=log2(2),
  y=3.2)
mtext("TMB",
  at=3,
  side=1,
  line=2.5,
  cex=axisCex)
mtext(paste("TMB", sig, sep=", "),
  at=3,
  side=3,
  line=1,
  cex=titleCex,
  font=titleF)

par(mar=oldMar)
addDots <- FALSE

sig <- "APOBEC3B"

oldMar <- par()$mar


## plot APOBEC expression by response
tmpDat <- pData(cds2)
tmpDat[, sig] <- m2[sig, ]
tmpDat <- tmpDat[tmpDat[, irf] != "NE",]
tmpDat[, irf] <- droplevels(tmpDat[, irf])

pval <- signif(wilcox.test(tmpDat[, sig] ~ tmpDat[, "binaryResponse"])$p.value, 
2)
print(paste("Wilcoxon test P APOBEC3B expression by binary response:", pval))
## [1] "Wilcoxon test P APOBEC3B expression by binary response: 0.0066"
pval <- signif(t.test(tmpDat[, sig] ~ tmpDat[, "binaryResponse"])$p.value, 
2)
print(paste("T test P APOBEC3B expression by binary response:", pval))
## [1] "T test P APOBEC3B expression by binary response: 0.0025"
yLab <- expression(paste("APOBEC3B expression (log"["2"], ")"))
par(mar=c(5.5, 4, 2, 0.5))

a <- boxplot(tmpDat[, sig] ~ tmpDat[, irf],
  ylab=yLab,
  col=color_palettes$irf_palette[levels(tmpDat[, irf])],
  cex.main=0.6,
  cex.axis=axisCex,
  cex.names=namesCex,
  cex.lab=labCex,  
  whisklty = 1)
if (addDots) {
  points(y=tmpDat[, sig],
    x=jitter(as.numeric(tmpDat[, irf]), factor=0.7),
    col=alpha("darkgrey", 0.7),
    pch=16)
}
mtext(a$n,
  side=3,
  at=1:nlevels(tmpDat[, irf]),
  line=0,
  cex=0.75)
mtext(paste(sig, "response", sep=", "),
  side=3,
  at=2.5,
  line=1,
  font=titleF,
  cex=titleCex)

## plot APOBEC expression by TMB
tmpDat <- pData(cds2)
tmpDat[, sig] <- m2[sig, ]
ind <- tmpDat[, ml] > 0
tmpDat[, ml] <- log2(tmpDat[, ml])

plot(tmpDat[ind, sig] ~ tmpDat[ind, ml],
    xlab="",
    ylab=expression(paste("APOBEC3B expression (log"["2"], ")")),
    pch=16,
    col="darkgrey",
    cex.axis=0.8,
    cex.lab=labCex,
    xaxt="n")
axis(1,
  at=log2(c(1,2,5,10,20,50)),
  labels=c(1,2,5,10,20,50),
  cex=axisCex)

fit <- lm(tmpDat[ind, sig] ~ tmpDat[ind, ml])
abline(fit, lty=2)

tmp <- cor.test(tmpDat[ind, sig], tmpDat[ind, ml], method="pearson")
text(paste("cor:", round(tmp$estimate, 2), 
    "\nP:", signif(tmp$p.value, 2)),
  font=2,
  cex=0.7,
  x=log2(2),
  y=2)
mtext("TMB",
  at=3,
  side=1,
  line=2.5,
  cex=axisCex)
mtext(paste("TMB", sig, sep=", "),
  at=3,
  side=3,
  line=1,
  cex=titleCex,
  font=titleF)

par(mar=oldMar)
path <- "CC Reg"

fmiPathP <- pathFMIout[pathFMIout$category == "TMB", ]
fmiGeneP <- pathFMIout3[pathFMIout3$category == "TMB", ]

ids <- ind_genes[[path]]
ids <- ids[ids %in% featureNames(fmi)]
fmiGoi <- fmi[featureNames(fmi) %in% ids, ]
sOrder <- order(as.numeric(pData(fmiGoi)[, ml]), decreasing=TRUE)

mutGoi <- any_mutation(fmiGoi)

tmp1 <- mutGoi[rowSums(mutGoi) > 0,]
tmp2 <- colSums(tmp1)
tmp2 <- tmp2 > 0
if ("TP53" %in% ids) {
  tmp3 <- mutGoi[ids[! ids %in% "TP53"], ]
  tmp3 <- colSums(tmp3)
  tmp3 <- tmp3 > 0
  tmp2 <- rbind(tmp2, tmp3)
}
toSort <- rowSums(tmp1)
toSort <- sort(toSort, decreasing=TRUE)
mutGoi <- rbind(tmp1[names(toSort),], tmp2)
rownames(mutGoi)[rownames(mutGoi) == "tmp2"] <- "pathway"
if ("TP53" %in% ids) {
  rownames(mutGoi)[rownames(mutGoi) == "tmp3"] <- "pathway w/o TP53"
}


mutGoi <- apply(mutGoi, 2, function(x) {
  ifelse(x, "mutant", "")
})

fmiCols <- c(mutant="black")
alter_fun = function(x, y, w, h, v) {
    if(v["mutant"]) grid.rect(x, y, w*0.9, h*0.9, gp = gpar(fill = fmiCols["mutant"], 
      col = NA))
  }


ha = HeatmapAnnotation(TMB = anno_barplot(as.numeric(pData(fmiGoi)[sOrder, ml]),
    border=FALSE,
    gp = gpar(fill="black")),
  annotation_legend_param=list(gp = gpar(fontsize = 1)),
  show_annotation_name = TRUE,
  annotation_name_side="left",
  annotation_name_gp = gpar(fontsize = 10))

rownames(mutGoi) <- paste0(rownames(mutGoi),
  ifelse(rownames(mutGoi) %in% fmiGeneP$Gene[fmiGeneP$AdjP < 0.05], "*", ""))

op1 <- oncoPrint(mutGoi[, sOrder], 
  col=fmiCols,
  alter_fun=alter_fun,
  row_order=NULL,
  column_order=NULL,
  column_title = "Mutations in cell cycle regulation genes", 
  top_annotation=ha,
  show_pct=TRUE,
  pct_gp = gpar(fontsize = 9),
  row_names_gp = gpar(fontsize = 9),
  column_title_gp = gpar(fontsize = 10),
  heatmap_legend_param=list(title ="Mutation Status"),
  show_row_barplot = FALSE,
  show_heatmap_legend = FALSE)
draw(op1, 
  padding = unit(c(1, 5, 0, 1), "mm"),
  row_title="mutation prevalence",
  row_title_gp=gpar(fontsize = 10))

ind <- !is.na(pData(fmi)$binaryResponse)
pData(fmi)$binaryResponse <- factor(pData(fmi)$binaryResponse,
  levels=c("CR/PR", "SD/PD"))

par(mar=c(5.5, 4.1, 2, 0))

path <- "CC Reg"

ids <- ind_genes[[path]]
ids <- ids[!ids %in% "TP53"]
tmp <- fmi[featureNames(fmi) %in% ids, ]
mutStatus <- any_mutation(tmp)
mutStatus <- colSums(mutStatus)

pval <- signif(fisher.test(table(pData(fmi)[ind, "binaryResponse"], mutStatus[ind] > 0))$p.value)
print(paste("Fisher P binary response by CCR mutation status:", pval))
## [1] "Fisher P binary response by CCR mutation status: 0.31104"
b <- table(responder = fmi$binaryResponse[ind],
        mutant = factor(mutStatus[ind] > 0))
nSamples <- colSums(b)
d <- prop.table(b, 
margin=2)

a <- barplot(d,
  legend.text=colnames(d),
  col=color_palettes$response_palette[rownames(d)],
  cex.names=namesCex,
  cex.axis=axisCex,
  cex.lab=labCex,
  width=0.2,
  xlim=c(0,1),
  ylab="fraction of patients",
  xaxt="n",
  args.legend=list(x="topright",
    bty="n",
    fill=color_palettes$response_palette,
    legend=c("CR/PR", "SD/PD"),
    cex=legendCex))
text(x = a, 
  y = par("usr")[3] - 0.06, 
  labels = ifelse(colnames(d) == FALSE, "non-mutant", "mutant"), 
  srt = -45, 
  xpd = TRUE,
  adj=0,
  cex=namesCex,
  xpd=TRUE)
mtext(nSamples,
  side=3,
  at=a,
  line=0,
  cex=0.7)
mtext("CCR w/o TP53, response",
  side=3,
  at=mean(a),
  line=1,
  font=titleF,
  cex=titleCex)

par(mar=oldMar)

ED Figure 2

orgMar <- par()$mar

keggRes <- read.csv(file.path(tablesDir,
    "SupplementaryTableS6.csv"),
  stringsAsFactor=FALSE,
  check.names=FALSE)
keggRes <- keggRes[, 1:6]
keggRes <- keggRes[!is.na(keggRes$S), ]
ind <- rev(1:nrow(keggRes))
keggRes <- keggRes[ind,]

up <- keggRes$Direction == "Up"
keggRes$"P (adj.)"[up] <- log10(keggRes$"P (adj.)"[up]) * -1
down <- keggRes$Direction == "Down"
keggRes$"P (adj.)"[down] <- log10(keggRes$"P (adj.)"[down])

resC1 <- read.csv(file.path(tablesDir,
    "SupplementaryTableS3.csv"),
  stringsAsFactor=FALSE,
  check.names=FALSE)

keggRes$RepGenes <- NA_character_
for (i in 1:nrow(keggRes)) {
  intGenes <- unlist(strsplit(keggRes$"Identified Genes"[i], ", "))
  if (length(intGenes) > 7) {
    intGenes <- resC1$Symbol[resC1$Symbol %in% intGenes][1:7]
    intGenes <- c(intGenes, "...")
  }
  intGenes <- paste(intGenes, collapse=", ")
  keggRes$RepGenes[i] <- intGenes
}

par(lheight=.7)
layout(matrix(c(1,0,0,0,1,0,0,0), 2, 4, byrow=TRUE))


proliferationPathways <- c("DNA replication", 
  "Cell cycle", 
  "Systemic lupus erythematosus", 
  "Alcoholism",
  "Oocyte meiosis",
  "Viral carcinogenesis",
  "Progesterone-mediated oocyte maturation")

setCols <- ifelse(keggRes$Name %in% proliferationPathways,
  "turquoise4",
    ifelse(keggRes$Name %in% c("Fanconi anemia pathway", "Homologous recombination", "Mismatch repair", "Nucleotide excision repair", "Base excision repair", "p53 signaling pathway"),
      "magenta4", 
        ifelse(keggRes$Name %in% "Cytokine-cytokine receptor interaction", "orange2", "darkgrey")))

a <- barplot(keggRes$"P (adj.)",
  horiz=TRUE,
  xlab="",
  xlim=c(-1, 17),
  col=setCols,
  axes=FALSE)
mtext("KEGG pathway",
  side=4,
  at=max(a[,1]) + 1.5,
  las=2,
  line=3,
  cex=titleCex)
mtext("Top-scoring genes",
  side=4,
  at=max(a[,1]) + 1.5,
  las=2,
  line=23,
  cex=titleCex)
mtext(side=4,
  at=a[,1],
  text=paste0(keggRes$Name,
    ifelse(keggRes$rep, "*", "")),
  las=2,
  line=1.5,
  cex=(namesCex - 0.1),
  col="black")
mtext(keggRes$RepGenes,
  side=4,
  at=a[,1],
  line=22,
  cex=(namesCex - 0.1),
  las=2)
axis(side=1,
  at=c(-2, seq(from=2, to=16, by=2)),
  labels=c(2, seq(from=2, to=16, by=2)),
  tick=TRUE,
  xpd=TRUE,
  cex=axisCex)
mtext(c("Up in responders", "Down"), 
  side=2, 
  line=1.5, 
  at=c(a[,1][10], a[,1][1]),
  cex=axisCex)
segments(x0=-2.1,
  x1=-2.1,
  y0=a[,1][2] - 0.3,
  y1=a[,1][length(a[,1])] + 0.3,
  col="black",
  xpd=TRUE,
  lwd=1)
segments(x0=-2.1,
  x1=-2.1,
  y0=a[,1][1] - 0.5,
  y1=a[,1][1] + 0.5,
  col="black",
  xpd=TRUE,
  lwd=1)

mtext(expression(paste("-log"["10"], " adj. p")),
  side=1,
  line=2.75,
  at=7,
  cex=axisCex)
legend("bottomright",
  fill=c("turquoise4", "magenta4", "orange2"),
  legend=c("Proliferation", "DNA damage response", expression(paste("TGF-", beta))),
  bty="n",
  cex=legendCex + 0.3)
segments(x0=0,
  x1=0,
  y0=-2,
  y1=max(a[,1])+1,
  col="black",
  xpd=TRUE,
  lwd=1)

par(mar=orgMar)
plotDots <- FALSE

orgMar <- par()$mar

tmpDat <- pData(cds2)
tmpDat$IFNG <- scale(m2["IFNG", ], center=TRUE, scale=TRUE)
tmpDat$IFNGR1 <- scale(m2["IFNGR1", ], center=TRUE, scale=TRUE)

tmpDat <- tmpDat[tmpDat[, irf] != "NE",]
tmpDat[, irf] <- droplevels(tmpDat[, irf])

for (sig in c("IFNG", "IFNGR1")) {

  yLim1 <- -3
  yLim2 <- ifelse(sig == "IFNG", 3, 5)
  yLim <- c(yLim1, yLim2)
  
  par(mar=c(5.5, 4.1, 2, 0.5))

  pval <- signif(wilcox.test(tmpDat[, sig] ~ tmpDat[, "binaryResponse"])$p.value, 
    2)
  print(paste("Wilcoxon test P", sig, "score by binary Response:", pval))

  pval <- signif(t.test(tmpDat[, sig] ~ tmpDat[, "binaryResponse"])$p.value, 
    2)
  print(paste("T test P", sig, "score by binary Response:", pval))

  yLab <- paste(sig, "expression")

  a <- boxplot(tmpDat[, sig] ~ tmpDat[, irf],
    ylab=yLab,
    #xlab="response",
    col=color_palettes[["irf_palette"]][levels(tmpDat[, irf])],
    cex.axis=axisCex,
    cex.names=namesCex,
    cex.lab=labCex, 
    ylim=yLim,
    whisklty = 1)
  if (plotDots) {
    points(y=tmpDat[, sig],
      x=jitter(as.numeric(tmpDat[, irf]), factor=0.7),
      col=alpha("darkgrey", 0.7),
      pch=16)
  }
  mtext(a$n,
    side=3,
    at=1:nlevels(tmpDat[, irf]),
    line=0,
    cex=0.75)
  mtext(paste0(sig, ", response"),
    side=3,
    at=2.5,
    line=1,
    font=titleF,
    cex=titleCex)

    yrange <- par("usr")[4] - par("usr")[3]
    yunit <- yrange/60
    segments(x0=1,
      x1=2,
      y0=par("usr")[4] - yunit * 4,
      y1=par("usr")[4] - yunit * 4,
      col="black",
      xpd=TRUE,
      lwd=1)
    segments(x0=3,
      x1=4,
      y0=par("usr")[4] - yunit * 4,
      y1=par("usr")[4] - yunit * 4,
      col="black",
      xpd=TRUE,
      lwd=1)
    segments(x0=1.5,
      x1=1.5,
      y0=par("usr")[4] - yunit * 4,
      y1=par("usr")[4] - yunit * 2,
      col="black",
      xpd=TRUE,
      lwd=1)
    segments(x0=3.5,
      x1=3.5,
      y0=par("usr")[4] - yunit * 4,
      y1=par("usr")[4] - yunit * 2,
      col="black",
      xpd=TRUE,
      lwd=1)
    segments(x0=1.5,
      x1=3.5,
      y0=par("usr")[4] - yunit * 2,
      y1=par("usr")[4] - yunit * 2,
      col="black",
      xpd=TRUE,
      lwd=1)
    text(y=par("usr")[4] - yunit ,
      x=2.5,
      labels=pLevel(pval),
      cex=ifelse(pLevel(pval) == "n.s.", axisCex - 0.1, axisCex)) 

}
## [1] "Wilcoxon test P IFNG score by binary Response: 5.2e-05"
## [1] "T test P IFNG score by binary Response: 9.1e-05"

## [1] "Wilcoxon test P IFNGR1 score by binary Response: 8.7e-05"
## [1] "T test P IFNGR1 score by binary Response: 0.00012"

par(mar=orgMar)
plotDots <- FALSE
binary <- FALSE

orgMar <- par()$mar

tmpDat <- pData(cds2)
tmpDat$TGFBR2 <- scale(m2["TGFBR2", ], center=TRUE, scale=TRUE)

tmpDat <- tmpDat[tmpDat[, irf] != "NE",]
tmpDat[, irf] <- droplevels(tmpDat[, irf])

sig <- "TGFBR2"
  

par(mar=c(5.5, 4.1, 2, 0.5))

pval <- signif(wilcox.test(tmpDat[, sig] ~ tmpDat[, "binaryResponse"])$p.value, 
  2)
print(paste("Wilcoxon test P", sig, "score by binary Response:", pval))
## [1] "Wilcoxon test P TGFBR2 score by binary Response: 0.00029"
pval <- signif(t.test(tmpDat[, sig] ~ tmpDat[, "binaryResponse"])$p.value, 
  2)
print(paste("T test P", sig, "score by binary Response:", pval))
## [1] "T test P TGFBR2 score by binary Response: 0.00019"
yLab <- paste(sig, "expression")

a <- boxplot(tmpDat[, sig] ~ tmpDat[, irf],
  ylab=yLab,
  col=color_palettes[["irf_palette"]][levels(tmpDat[, irf])],
  cex.axis=axisCex,
  cex.names=namesCex,
  cex.lab=labCex, 
  ylim=c(-5, 5),
  whisklty = 1)
if (plotDots) {
  points(y=tmpDat[, sig],
    x=jitter(as.numeric(tmpDat[, irf]), factor=0.7),
    col=alpha("darkgrey", 0.7),
    pch=16)
}
mtext(a$n,
  side=3,
  at=1:nlevels(tmpDat[, irf]),
  line=0,
  cex=0.75)
mtext(paste0(sig, ", response"),
  side=3,
  at=2.5,
  line=1,
  font=titleF,
  cex=titleCex)

yrange <- par("usr")[4] - par("usr")[3]
yunit <- yrange/60
segments(x0=1,
  x1=2,
  y0=par("usr")[4] - yunit * 4,
  y1=par("usr")[4] - yunit * 4,
  col="black",
  xpd=TRUE,
  lwd=1)
segments(x0=3,
  x1=4,
  y0=par("usr")[4] - yunit * 4,
  y1=par("usr")[4] - yunit * 4,
  col="black",
  xpd=TRUE,
  lwd=1)
segments(x0=1.5,
  x1=1.5,
  y0=par("usr")[4] - yunit * 4,
  y1=par("usr")[4] - yunit * 2,
  col="black",
  xpd=TRUE,
  lwd=1)
segments(x0=3.5,
  x1=3.5,
  y0=par("usr")[4] - yunit * 4,
  y1=par("usr")[4] - yunit * 2,
  col="black",
  xpd=TRUE,
  lwd=1)
segments(x0=1.5,
  x1=3.5,
  y0=par("usr")[4] - yunit * 2,
  y1=par("usr")[4] - yunit * 2,
  col="black",
  xpd=TRUE,
  lwd=1)
text(y=par("usr")[4] - yunit ,
  x=2.5,
  labels=pLevel(pval),
  cex=ifelse(pLevel(pval) == "n.s.", axisCex - 0.1, axisCex))   

par(mar=orgMar)
bio <- "TGFBR2"

tmpDat  <- pData(cds2)
tmpDat$TGFBR2 <- m2["TGFBR2", ]

tmpDat <- tmpDat[!is.na(tmpDat[, bio]),]
tmpDat$group <- cut(tmpDat[, bio], 
  quantile(tmpDat[, bio], c(0, 0.25, 0.5, 0.75, 1)),
  include.lowest = TRUE)
tmpDat$group <- factor(tmpDat$group,
  labels=c("1st Quart", 
    "2nd Quart",
    "3rd Quart",
    "4th Quart"))
fittedSurv <- survfit(Surv(time=os, event=censOS) ~ group, 
              data=tmpDat)
diffSurv <- survdiff(Surv(time=os, event=censOS) ~ group,           
            data=tmpDat)
plotSurvival2(survFit=fittedSurv, 
    survDiff=diffSurv, 
    diff.factor=as.factor(tmpDat$group), 
    main="TGFBR2, survival", 
    cols=rev(color_palettes[["irf_palette"]][1:4]),
    ltypes=rep(1, nlevels(tmpDat$group)),
    pval="none",
    #coxphData=coxphDat,
    plotMedian=TRUE,
    plot.nrisk=FALSE,
    xLab="OS (months)",
    lwd=3,
    cexMedian=axisCex - 0.1,
    cexLegend=legendCex,
    cex.lab=labCex)

print(coxph(Surv(os, censOS) ~ 
  as.integer(tmpDat$group),
  data=tmpDat))
## Call:
## coxph(formula = Surv(os, censOS) ~ as.integer(tmpDat$group), 
##     data = tmpDat)
## 
##                             coef exp(coef) se(coef)     z      p
## as.integer(tmpDat$group) 0.13451   1.14398  0.05866 2.293 0.0218
## 
## Likelihood ratio test=5.29  on 1 df, p=0.02147
## n= 348, number of events= 232
goi <- c("CD 8 T effector",
  "gene19")

resp <- "binaryResponse"
tmpDat <- pData(cds2)[!is.na(pData(cds2)[, resp]), ]
tmpDat[, resp] <- droplevels(tmpDat[, resp])
for (sig in goi) {
  tmpDat[, sig] <- scale(tmpDat[, sig], center=TRUE, scale=TRUE) 
}

tmpDat <- tmpDat[!is.na(tmpDat[, ml]) & tmpDat[, ml] > 0, ]
tmpDat[, ml] <- log2(tmpDat[, ml])

tmpDat <- tmpDat[! is.na(tmpDat$"Immune phenotype"),]

tmpDat$isExcluded <- as.factor(ifelse(tmpDat[, "Immune phenotype"] == "excluded",
  "Yes", "No"))

varExp1 <- t(sapply(c(goi, ml), function(sig) {
  fit <- glm(tmpDat[, resp] ~ 
    tmpDat[, sig],
    family="binomial")
  fit$null.deviance - fit$deviance
}))
varExp1 <- setNames(varExp1[1,],
  colnames(varExp1))

fit2.2 <- glm(tmpDat[, resp] ~ 
    tmpDat[, "CD 8 T effector"] + tmpDat[, "gene19"],
    family="binomial")
varExp2.2 <- fit2.2$null.deviance - fit2.2$deviance

fit2.4 <- glm(tmpDat[, resp] ~ 
    tmpDat[, ml] + tmpDat[, "CD 8 T effector"], 
    family="binomial")
varExp2.4 <- fit2.4$null.deviance - fit2.4$deviance

fit2.6 <- glm(tmpDat[, resp] ~ 
    tmpDat[, ml] * tmpDat[, "gene19"], 
    family="binomial")
varExp2.6 <- fit2.6$null.deviance - fit2.6$deviance

varExp2 <- c(varExp2.2, varExp2.4, varExp2.6)
names(varExp2) <- c("Teff,TBRS", "Teff,TMB", "TBRS,TMB")
varExp2 <- sort(varExp2)

fit3.2 <- glm(tmpDat[, resp] ~ 
    tmpDat[, "CD 8 T effector"] + tmpDat[, ml] + tmpDat[, "gene19"] +
    tmpDat[, ml]:tmpDat[, "gene19"] +
    tmpDat[, "gene19"]*tmpDat$isExcluded,
    family="binomial")
varExp3.2 <- fit3.2$null.deviance - fit3.2$deviance
names(varExp3.2) <- "Teff,TBRS,TMB"

allVars <- c(varExp1, 
    NA, 
    varExp2, 
    NA, 
    varExp3.2)
names(allVars) <- sub("CD 8 T effector", "Teff", names(allVars))
names(allVars) <- sub("gene19", "TBRS", names(allVars))
names(allVars) <- sub(ml, "TMB", names(allVars))

## repeat previous model to fit displayed data
fit_teff_0 <- glm(tmpDat[, resp] ~ 
  tmpDat[, "CD 8 T effector"],
  family="binomial")
fit_tgfb_0 <- glm(tmpDat[, resp] ~ 
  tmpDat[, "gene19"],
  family="binomial")
fit_teff_tgfb <- glm(tmpDat[, resp] ~ 
  tmpDat[, "CD 8 T effector"] +
  tmpDat[, "gene19"],
  family="binomial")

fit_tt_0 <- glm(tmpDat[, resp] ~ 
  tmpDat[, "CD 8 T effector"] +
  tmpDat[, "gene19"],
  family="binomial")
fit_ml_0 <- glm(tmpDat[, resp] ~ 
  tmpDat[, ml],
  family="binomial")
fit_tt_ml <- glm(tmpDat[, resp] ~ 
  tmpDat[, "CD 8 T effector"] +
  tmpDat[, "gene19"] +
  tmpDat[, ml],
  family="binomial")
fit_tgfb_ml <- glm(tmpDat[, resp] ~ 
  tmpDat[, "gene19"] +
  tmpDat[, ml],
  family="binomial")
fit_tgfb_ml_int <- glm(tmpDat[, resp] ~ 
  tmpDat[, "gene19"] *
  tmpDat[, ml],
  family="binomial")
fit_teff_ml <- glm(tmpDat[, resp] ~ 
  tmpDat[, "CD 8 T effector"] +
  tmpDat[, ml],
  family="binomial")
fit_core_int <- glm(tmpDat[, resp] ~ 
  tmpDat[, "CD 8 T effector"] +
  tmpDat[, ml] +
  tmpDat[, "gene19"] +
  tmpDat[, "gene19"]:tmpDat[, ml] +
  tmpDat[, "gene19"]*tmpDat$isExcluded,
  family="binomial")

# likelihood ratio test Ps for display
tgfbP <- signif(anova(fit_tgfb_0, fit_teff_tgfb, test="Chisq")$"Pr(>Chi)"[2], 2)
teffP <- signif(anova(fit_teff_0, fit_teff_tgfb, test="Chisq")$"Pr(>Chi)"[2], 2)
ttP <- signif(anova(fit_tt_0, fit_tt_ml, test="Chisq")$"Pr(>Chi)"[2], 2)
mlP <- signif(anova(fit_ml_0, fit_tt_ml, test="Chisq")$"Pr(>Chi)"[2], 2)

teffmlP <- signif(anova(fit_ml_0, fit_teff_ml, test="Chisq")$"Pr(>Chi)"[2], 2)
teffmlP2 <- signif(anova(fit_teff_0, fit_teff_ml, test="Chisq")$"Pr(>Chi)"[2], 2)

tgfbmlP <- signif(anova(fit_ml_0, fit_tgfb_ml, test="Chisq")$"Pr(>Chi)"[2], 2)
tgfbmlP2 <- signif(anova(fit_tgfb_0, fit_tgfb_ml, test="Chisq")$"Pr(>Chi)"[2], 2)
tgfbmlP3 <- signif(anova(fit_tgfb_ml, fit_tgfb_ml_int, test="Chisq")$"Pr(>Chi)"[2], 2)

tgfbmlP_int <- signif(anova(fit_tgfb_0, fit_tgfb_ml_int, test="Chisq")$"Pr(>Chi)"[2], 2)
tgfbmlP_int2 <- signif(anova(fit_ml_0, fit_tgfb_ml_int, test="Chisq")$"Pr(>Chi)"[2], 2)

tt_tgfb_mlP <- signif(anova(fit_tt_ml, fit_tgfb_ml, test="Chisq")$"Pr(>Chi)"[2], 2)

core_interaction <- signif(anova(fit_tt_ml, fit_core_int, test="Chisq")$"Pr(>Chi)"[2], 2)

core_tt <- signif(anova(fit_tt_0, fit_core_int, test="Chisq")$"Pr(>Chi)"[2], 2)
core_tgfbml <- signif(anova(fit_tgfb_ml_int, fit_core_int, test="Chisq")$"Pr(>Chi)"[2], 2)
core_teffml <- signif(anova(fit_teff_ml, fit_core_int, test="Chisq")$"Pr(>Chi)"[2], 2)

pvals <- list(tgfbP=tgfbP,
  teffP=teffP,
  ttP=ttP,
  mlP=mlP,
  teffmlP=teffmlP,
  teffmlP2=teffmlP2,
  tgfbmlP=tgfbmlP,
  tgfbmlP2=tgfbmlP2,
  tgfbmlP3=tgfbmlP3,
  tt_tgfb_mlP=tt_tgfb_mlP,
  core_int=core_interaction,
  tgfbmlP_int=tgfbmlP_int,
  tgfbmlP_int2=tgfbmlP_int2,
  core_tt=core_tt,
  core_tgfbml=core_tgfbml,
  core_teffml=core_teffml)
print(pvals)
## $tgfbP
## [1] 0.0032
## 
## $teffP
## [1] 0.0026
## 
## $ttP
## [1] 9e-07
## 
## $mlP
## [1] 0.08
## 
## $teffmlP
## [1] 0.2
## 
## $teffmlP2
## [1] 4.9e-08
## 
## $tgfbmlP
## [1] 0.23
## 
## $tgfbmlP2
## [1] 6.6e-08
## 
## $tgfbmlP3
## [1] 0.014
## 
## $tt_tgfb_mlP
## [1] 0.057
## 
## $core_int
## [1] 0.0052
## 
## $tgfbmlP_int
## [1] 2.2e-08
## 
## $tgfbmlP_int2
## [1] 0.024
## 
## $core_tt
## [1] 1.9e-07
## 
## $core_tgfbml
## [1] 0.016
## 
## $core_teffml
## [1] 0.0028
pvals <- lapply(pvals, function(x) {
  ifelse(x < 0.001, "***",
  ifelse(x < 0.01, "**",
    ifelse(x < 0.05, "*",
      ifelse(x < 0.1, ".", "n.s."))))
})


oldMar <-  par()$mar

par(mar=c(6, 3.9, 2, 0))

a <- barplot(allVars,
  main="",
  ylab="explained variance (%)",
  col="grey",
  width=0.06,
  xlim=c(0,1),
  ylim=c(0,58),
  #args.legend=list(bg=alpha("white", 0.7),
  #  box.col=alpha("white", 0.7),
  #  cex=0.9),
  xaxt="n",
  cex.axis=0.9)
axis(1, 
  at=a[c(1:3, 5:7, 9)],
  labels=FALSE)
text(x = a, 
  y=-3, 
  labels = names(allVars),
  srt = -45, 
  xpd = TRUE,
  adj=0,
  cex=namesCex)
print("variance explained:")
## [1] "variance explained:"
print(allVars)
##          Teff          TBRS           TMB                   Teff,TBRS      Teff,TMB 
##      4.125469      4.491943     32.246450            NA     13.169119     33.885674 
##      TBRS,TMB               Teff,TBRS,TMB 
##     39.742273            NA     50.042237
mtext(at = mean(a), 
  line=0.85, 
  side=3,
  text="Variance explained, core pathways",
  cex=titleCex,
  font=titleF)

## improvement of Teff, TBRS over single models
text(x = a[5], 
  y=c(16.3, 14.8), 
  labels = unlist(pvals[c("teffP", "tgfbP")]), 
  #pos=3,
  cex=0.9,
  xpd=TRUE)

## improvement of Teff, TMB over single models
text(x = a[6], 
  y=c(37.3, 35.5), 
  labels = unlist(pvals[c("teffmlP2", "teffmlP")]), 
  #pos=3,
  cex=0.9,
  xpd=TRUE)

## improvement of TBRS, TMB over single models
text(x = a[7], 
  y=c(42.5, 41), 
  labels = unlist(pvals[c("tgfbmlP_int", "tgfbmlP_int2")]), 
  #pos=3,
  cex=0.9,
  xpd=TRUE)

## improvement of double on final core
segments(x0=a[7],
  x1=a[9],
  y0=51,
  y1=51,
  col="black",
  xpd=TRUE,
  lwd=2)
text(x = median(c(a[7], a[9])), 
  y=52.5, 
  labels = pvals[["core_tgfbml"]], 
  #pos=3,
  cex=0.9,
  xpd=TRUE)
segments(x0=a[6],
  x1=a[9],
  y0=55,
  y1=55,
  col="black",
  xpd=TRUE,
  lwd=2)
text(x = median(c(a[6], a[9])), 
  y=56.5, 
  labels = pvals[["core_teffml"]], 
  #pos=3,
  cex=0.9,
  xpd=TRUE)
segments(x0=a[5],
  x1=a[9],
  y0=59,
  y1=59,
  col="black",
  xpd=TRUE,
  lwd=2)
text(x = median(c(a[5], a[9])), 
  y=60.5, 
  labels = pvals[["core_tt"]], 
  #pos=3,
  cex=0.9,
  xpd=TRUE)

par(mar=oldMar)



Session Info

[1] “2019-11-16 23:37:15” R version 3.5.3 (2019-03-11) Platform: x86_64-apple-darwin15.6.0 (64-bit) Running under: macOS Mojave 10.14.4

Matrix products: default BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib LAPACK: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRlapack.dylib

locale: [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages: [1] stats4 parallel grid stats graphics grDevices utils datasets [9] methods base

other attached packages: [1] corrplot_0.84 DT_0.9
[3] DESeq2_1.22.2 SummarizedExperiment_1.12.0
[5] DelayedArray_0.8.0 BiocParallel_1.16.6
[7] matrixStats_0.55.0 GenomicRanges_1.34.0
[9] GenomeInfoDb_1.18.2 IRanges_2.16.0
[11] S4Vectors_0.20.1 rmarkdown_1.16
[13] knitr_1.25 biomaRt_2.38.0
[15] readxl_1.3.1 edgeR_3.24.3
[17] limma_3.38.3 ROCR_1.0-7
[19] gplots_3.0.1.1 broom_0.5.2
[21] glmnet_2.0-18 foreach_1.4.7
[23] Matrix_1.2-17 MASS_7.3-51.4
[25] caret_6.0-84 lattice_0.20-38
[27] IMvigor210CoreBiologies_1.0.0 Biobase_2.42.0
[29] BiocGenerics_0.28.0 survminer_0.4.6
[31] survival_3.1-7 seriation_1.2-8
[33] viridis_0.5.1 viridisLite_0.3.0
[35] ggbeeswarm_0.6.0 ggpubr_0.2.3
[37] ggsci_2.9 GenomicDataCommons_1.6.0
[39] magrittr_1.5 PerformanceAnalytics_1.5.3
[41] xts_0.11-2 zoo_1.8-6
[43] ComplexHeatmap_1.20.0 forcats_0.4.0
[45] stringr_1.4.0 dplyr_0.8.3
[47] purrr_0.3.3 readr_1.3.1
[49] tidyr_1.0.0 tibble_2.1.3
[51] ggplot2_3.2.1 tidyverse_1.2.1
[53] workflowr_1.5.0

loaded via a namespace (and not attached): [1] utf8_1.1.4 tidyselect_0.2.5 htmlwidgets_1.5.1
[4] RSQLite_2.1.2 AnnotationDbi_1.44.0 TSP_1.1-7
[7] DESeq_1.34.1 munsell_0.5.0 codetools_0.2-16
[10] withr_2.1.2 colorspace_1.4-1 rstudioapi_0.10
[13] ggsignif_0.6.0 labeling_0.3 git2r_0.26.1
[16] GenomeInfoDbData_1.2.0 KMsurv_0.1-5 bit64_0.9-7
[19] rprojroot_1.3-2 vctrs_0.2.0 generics_0.0.2
[22] ipred_0.9-9 xfun_0.10 gclus_1.3.2
[25] R6_2.4.0 locfit_1.5-9.1 bitops_1.0-6
[28] assertthat_0.2.1 promises_1.1.0 scales_1.0.0
[31] nnet_7.3-12 beeswarm_0.2.3 gtable_0.3.0
[34] processx_3.4.1 timeDate_3043.102 rlang_0.4.1
[37] zeallot_0.1.0 genefilter_1.64.0 GlobalOptions_0.1.1
[40] splines_3.5.3 lazyeval_0.2.2 acepack_1.4.1
[43] ModelMetrics_1.2.2 checkmate_1.9.4 reshape2_1.4.3
[46] BiocManager_1.30.9 yaml_2.2.0 modelr_0.1.5
[49] crosstalk_1.0.0 backports_1.1.5 httpuv_1.5.2
[52] Hmisc_4.3-0 lava_1.6.6 tools_3.5.3
[55] ellipsis_0.3.0 RColorBrewer_1.1-2 Rcpp_1.0.3
[58] plyr_1.8.4 base64enc_0.1-3 progress_1.2.2
[61] zlibbioc_1.28.0 RCurl_1.95-4.12 ps_1.3.0
[64] prettyunits_1.0.2 rpart_4.1-15 GetoptLong_0.1.7
[67] cowplot_1.0.0 haven_2.2.0 cluster_2.1.0
[70] fs_1.3.1 data.table_1.12.6 circlize_0.4.8
[73] mime_0.7 hms_0.5.2 evaluate_0.14
[76] xtable_1.8-4 XML_3.98-1.20 gridExtra_2.3
[79] shape_1.4.4 compiler_3.5.3 KernSmooth_2.23-16
[82] crayon_1.3.4 htmltools_0.4.0 later_1.0.0
[85] Formula_1.2-3 geneplotter_1.60.0 lubridate_1.7.4
[88] DBI_1.0.0 rappdirs_0.3.1 cli_1.1.0
[91] quadprog_1.5-7 gdata_2.18.0 gower_0.2.1
[94] pkgconfig_2.0.3 km.ci_0.5-2 registry_0.5-1
[97] foreign_0.8-72 recipes_0.1.7 xml2_1.2.2
[100] annotate_1.60.1 vipor_0.4.5 XVector_0.22.0
[103] prodlim_2019.10.13 rvest_0.3.5 callr_3.3.2
[106] digest_0.6.22 cellranger_1.1.0 htmlTable_1.13.2
[109] survMisc_0.5.5 dendextend_1.12.0 curl_4.2
[112] shiny_1.4.0 gtools_3.8.1 rjson_0.2.20
[115] lifecycle_0.1.0 nlme_3.1-142 jsonlite_1.6
[118] fansi_0.4.0 pillar_1.4.2 fastmap_1.0.1
[121] httr_1.4.1 glue_1.3.1 iterators_1.0.12
[124] bit_1.1-14 class_7.3-15 stringi_1.4.3
[127] blob_1.2.0 latticeExtra_0.6-28 caTools_1.17.1.2
[130] memoise_1.1.0